---
title: "Assessing Canada’s Drug Shortage Problem"
output: 
  flexdashboard::flex_dashboard:
    orientation: columns
    vertical_layout: fill
    source_code: embed
---

```{r setup, include=FALSE}
library(flexdashboard)
library(tidyverse, warn.conflicts = FALSE)
library(lubridate)

# settings
FAR_FUTURE = parse_datetime("2100-12-31")
FAR_PAST = parse_datetime("1970-01-01")
DATA_UPDATED_DATE = parse_datetime("2019-12-25")

DPD_DATA_DIR = 'data/dpd/'  # Load the Health Canada Drug Product Database
DPD_DOWNLOAD = FALSE   # set to TRUE to attempt to download the latest DPD dataset
DPD_URL_BASE='http://www.hc-sc.gc.ca/dhp-mps/alt_formats/zip/prodpharma/databasdon/'

theme_set(theme_minimal())
caption = list(labs(caption =  paste(
  paste("source: drugshortagescanada.ca, last updated:", DATA_UPDATED_DATE), 
  "jon.pipitone.ca/research/drug-shortages", 
  sep="\n")))

# load data
DSCDATA= "data/dsc_shortages_2020-12-29.csv"
dsc_df = read_csv(DSCDATA, trim_ws = T, 
                  col_types = cols(
                    status = col_factor(),
                    shortage_reason_en = col_factor()
                  )) %>%
  select(-starts_with("fr_"), -ends_with("_fr")) %>%
  mutate(
    reason = fct_recode(shortage_reason_en,  
      "Demand increase" = "Demand increase for the drug.", 
      "Manufacturing distruption" = "Disruption of the manufacture of the drug.",
      "Ingredient shortage" = "Shortage of an active ingredient.", 
      "Ingredient shortage" = "Shortage of an inactive ingredient or component.",
      "Manufacturing practices" = "Requirements related to complying with good manufacturing practices.",
      "Other" = "Other (Please describe in comments)",
      "Shipping delay" = "Delay in shipping of the drug." 
    ) 
  ) %>%
  mutate(
    end_date_no_na = coalesce(actual_end_date, estimated_end_date, FAR_FUTURE), 
    start_date_no_na = coalesce(actual_start_date, anticipated_start_date, FAR_PAST))

## Parsual EML
# https://dx.doi.org/10.9778%2Fcmajo.20160122
emlist = read_csv('data/cleanmeds.csv', trim_ws=T, col_types = cols()) %>% 
  select(category, medication, ATC)

## DPD
if(DPD_DOWNLOAD) {
  for (f in c('allfiles.zip', 'allfiles_ia.zip', 'allfiles_ap.zip', 'allfiles_dr.zip')) {
    destfile = paste0(DPD_DATA_DIR, f)
    download.file(paste0(DPD_URL_BASE, f), destfile)
    unzip(destfile, exdir = DPD_DATA_DIR)
  }
}

tables = list(
  comp = c("DRUG_CODE","MFR_CODE","COMPANY_CODE","COMPANY_NAME","COMPANY_TYPE",
           "ADDRESS_MAILING_FLAG","ADDRESS_BILLING_FLAG",
           "ADDRESS_NOTIFICATION_FLAG","ADDRESS_OTHER","SUITE_NUMBER",
           "STREET_NAME","CITY_NAME","PROVINCE","COUNTRY","POSTAL_CODE",
           "POST_OFFICE_BOX","PROVINCE_F","COUNTRY_F"),
  form = c("DRUG_CODE","PHARM_FORM_CODE","PHARMACEUTICAL_FORM",
           "PHARMACEUTICAL_FORM_F"),
  ingred = c("DRUG_CODE","ACTIVE_INGREDIENT_CODE","INGREDIENT",
             "INGREDIENT_SUPPLIED_IND","STRENGTH","STRENGTH_UNIT",
             "STRENGTH_TYPE","DOSAGE_VALUE","BASE","DOSAGE_UNIT","NOTES",
             "INGREDIENT_F","STRENGTH_UNIT_F","STRENGTH_TYPE_F","DOSAGE_UNIT_F"),
  package = c("DRUG_CODE","UPC","PACKAGE_SIZE_UNIT","PACKAGE_TYPE",
              "PACKAGE_SIZE","PRODUCT_INFORMATION","PACKAGE_SIZE_UNIT_F",
              "PACKAGE_TYPE_F"),
  pharm = c("DRUG_CODE","PHARMACEUTICAL_STD"),
  drug = c("DRUG_CODE","PRODUCT_CATEGORIZATION","CLASS",
           "DRUG_IDENTIFICATION_NUMBER","BRAND_NAME","DESCRIPTOR",
           "PEDIATRIC_FLAG","ACCESSION_NUMBER","NUMBER_OF_AIS",
           "LAST_UPDATE_DATE","AI_GROUP_NO","CLASS_F","BRAND_NAME_F",
           "DESCRIPTOR_F"),
  route = c("DRUG_CODE","ROUTE_OF_ADMINISTRATION_CODE",
            "ROUTE_OF_ADMINISTRATION","ROUTE_OF_ADMINISTRATION_F"),
  sched = c("DRUG_CODE","SCHEDULE","SCHEDULE_F"),
  status = c("DRUG_CODE","CURRENT_STATUS_FLAG","STATUS","HISTORY_DATE",
             "STATUS_F","LOT_NUMBER","EXPIRATION_DATE"),
  ther = c("DRUG_CODE","TC_ATC_NUMBER","TC_ATC","TC_AHFS_NUMBER","TC_AHFS",
           "TC_ATC_F","TC_AHFS_F"),
  vet = c("DRUG_CODE","VET_SPECIES","VET_SUB_SPECIES","VET_SPECIES_F")
)

for (table_name in names(tables)) {
  headers = tables[[table_name]]
  # create new variables based on the DPD table names (e.g. dpd_vet) and load in the
  # data from the corresponding files (active, inactive, etc..) and merge them
  assign(paste0('dpd_', table_name), 
    dir(DPD_DATA_DIR, pattern = paste0(table_name, '.*.txt$'), full.names = TRUE) %>%
      keep(function(x) file.size(x) > 0) %>%
      map(read_csv, 
          col_names = headers, 
          locale = locale('en', date_format = '%d-%b-%Y'), 
          col_types = map(headers,  function (x) # guess base on name 
            case_when(
              endsWith(x, '_DATE') ~ 'D', 
              x == 'DRUG_CODE' ~ 'i', 
              TRUE ~ 'c')) %>% 
          setNames(headers), progress = FALSE) %>%
      reduce(rbind) %>%
      select(-ends_with("_F"))  # je ne parle pas bien le francais 
  )
}

## Add the DPD ATC number to the DSC data
din_to_atc = inner_join(dpd_drug, dpd_ther, by="DRUG_CODE") %>% 
  select(din = DRUG_IDENTIFICATION_NUMBER, dpd_atc_number = TC_ATC_NUMBER) %>% 
  distinct()

dsc_df = left_join(dsc_df, din_to_atc, by=c("din"))
```

Shortages
==========

Column
-----------------------------------------------------------------------

### Total number of drugs in shortage over time
```{r}
time.by.month = tibble(
  month = seq(parse_date("2017-01-01"), parse_date("2020-01-01"), by = "1 month"), 
  i = 1
)
# https://community.rstudio.com/t/tidy-way-to-range-join-tables-on-an-interval-of-dates/7881
shortages.by.month = 
  dsc_df %>%
  mutate(i = 1) %>%
  full_join(time.by.month, by = "i") %>%
  filter(int_overlaps(
    interval(month, month+months(1)-days(1)),
    interval(floor_date(actual_start_date, "month"), end_date_no_na))) %>%
  group_by(month) %>%
  summarize(n = n())

ggplot(shortages.by.month, aes(x = month, y = n)) + 
  geom_point() + geom_line() + 
  labs(x = "Date", y = "Number of drugs in shortage") + 
  caption
```

### Shortage start date

```{r, fig.asp=0.70}
dsc_df %>% 
  mutate(Month = floor_date(actual_start_date, "month")) %>% 
  filter(between(actual_start_date, parse_datetime("2017-01-01"), parse_datetime("2020-01-01")), 
         status != "avoided_shortage") %>%
  group_by(Month, status) %>% 
  summarize(count = n()) %>%
  ggplot(aes(x=Month, y = count, color = status)) + 
    geom_point() + geom_line() +  
    #geom_col(position="stack") + 
    scale_x_datetime(
      date_breaks = "1 year", 
      date_labels = "%Y") + 
    theme(legend.position="bottom") +
    labs(x = "Shortage start date", y = "Number of shortages", fill = "Current status") + 
    caption
```

Column 
-----------------------------------------------------------------------

### Duration of shortages
```{r}
data = dsc_df %>%
  filter(status %in% c("active_confirmed", "resolved")) %>%
  mutate(duration = interval(actual_start_date, actual_end_date)/weeks(1), 
         year = year(actual_start_date)) %>%
  filter(duration > 0, year >= 2017)

data %>%
  ggplot(aes(x = ceiling(duration))) + 
    geom_histogram(bins=50, position="identity") + 
    facet_grid(. ~ year) + 
    labs(fill = "Year", x = "Duration (weeks)", y = "# of shortages") + 
    theme(legend.position = "none") + 
    caption
```


### Causes of drug shortages

```{r, fig.asp=0.70}
dsc_df %>% 
  filter(status %in% c("resolved", "active_confirmed")) %>%
  mutate(reason=fct_infreq(reason)) %>%
  group_by(reason) %>%
  summarize(shortages=n()) %>%
  mutate(percent=shortages/sum(shortages)*100) %>% 
  ggplot(aes(x=2,y=percent,fill=reason)) + 
  xlim(0.5,2.5) + 
  geom_bar(width=1,stat="identity") + 
  coord_polar("y") + 
  geom_text(aes(label = paste0(round(percent), "%")), position = position_stack(vjust = 0.5), size=3) +
  labs(x=NULL,y=NULL,fill=NULL,
       title="Causes of drug shortages") + 
  theme_classic() + 
  theme(
    axis.line=element_blank(), 
    axis.text=element_blank(), 
    axis.ticks=element_blank()) +
  caption
```

Column {data-width=350}
-----------------------------------------------------------------------

### Active shortages
```{r}
active_shortages = dsc_df %>% 
  filter(
    status == "active_confirmed", 
    now() %within% interval(actual_start_date, end_date_no_na))
valueBox(nrow(active_shortages), icon="fa-capsules")
```

### Anticipated shortages next week
```{r}
week_start = floor_date(now(), "week")
this_week = interval(week_start, week_start + weeks(1))
next_week = int_shift(this_week, by = weeks(1))
expected_shortages = dsc_df %>% 
  filter(
    int_overlaps(next_week, interval(actual_start_date, end_date_no_na)))
valueBox(nrow(expected_shortages), icon="fa-clock")
```

### Summary (via the database)

```{r}
knitr::kable(dsc_df %>% group_by(status) %>% summarize(n = n()))
```


Essential Medicines 
=======================

Column 
-----------------

### Essential medicines impacted by shortages in 2019
```{r}
emls_with_shortages = dsc_df %>%
  filter(int_overlaps("2019-01-01" %--% "2019-12-31", actual_start_date %--% end_date_no_na)) %>%
  right_join(emlist, by=c("dpd_atc_number" = "ATC"))  %>%
  filter(!is.na(din))  %>%
  group_by(category, medication, dpd_atc_number) %>%
  summarize(Shortages = n()) %>%
  rename(Category = category, 
         Medication = medication,
         ATC = dpd_atc_number) %>%
  arrange(desc(Shortages))
knitr::kable(emls_with_shortages)
```

### All 2019 shortages by EML category
```{r}
shortages_by_eml = dsc_df %>%
  filter(int_overlaps("2019-01-01" %--% "2019-12-31", actual_start_date %--% end_date_no_na)) %>%
  left_join(emlist, by=c("dpd_atc_number" = "ATC")) %>%
  mutate(category = replace_na(category, "Not an essential medicine")) %>%
  mutate(category = str_replace(category, " -.*", "")) %>%
  group_by(category) %>%
  summarize(Shortages = n()) %>%
  mutate(percent_shortages = round(Shortages/sum(Shortages) * 100, digits=2)) %>%
  arrange(desc(percent_shortages)) %>%
  rename(Category = category, 
         "% of all shortages" = percent_shortages)
knitr::kable(shortages_by_eml)
```

Column {data-width=250}
-----------------------------------------------------------------------

### Shortages of Essential Medicines in 2019
```{r}
valueBox(nrow(emls_with_shortages), icon="fa-ambulance")
```

### Of all shortages in 2019
```{r}
valueBox(
  paste0(100-shortages_by_eml$'% of all shortages'[[1]], "%"),
  icon="fa-chart-pie")
```





About
=====

Column
------

### About this website

Hello. I am Jon Pipitone, medical student at Queen's University, Kingston, Ontario, Canada.

This website is a work-in-progress and is a place to display my ongoing research into drug shortages in Canada. You can follow my progress on my blog https://jon.pipitone.ca/blog. 

This work is done under supervision by [Dr. Jacalyn Duffin](https://jacalynduffin.ca/), and is an extension of our earlier analysis: 

*Donelle, Jessy, Jacalyn Duffin, Jonathan Pipitone, and Brian White-Guay. “Assessing Canada’s Drug Shortage Problem.” CD Howe Institute Commentary 515, 2018. https://doi.org/10.2139/ssrn.3192558.*

All of the code that powers this site, and other goodies can be found at https://github.com/pipitone/drug-shortages

*Disclaimer:* As this is ongoing research, the analysis and results presented here is likely to change, and has not undergone any peer review. 

Questions or comments? Please send me an email at jon@pipitone.ca

Column  {data-width=100}
----------------------------

### Data sources
##### [drugshortagescanada.ca](https://www.drugshortagescanada.ca)
Database of shortage and discontinuation reports.  Active from March 2017 to
present. Manufacturers are required to report.

##### [Drug Product Database](https://www.canada.ca/en/health-canada/services/drugs-health-products/drug-products/drug-product-database.html)
Database of historical records of drug products in Canada. Maintained by Health Canada.

##### [CLEANMeds Essential Medicines List](https://cleanmeds.ca)

### Terminology
**DSC**: Drug Shortages Canada database

**DPD**: Health Canada's Drug Products Database

**ATC**: Anatomical Therapeutic Chemical Classification system